Brownian motion and Diffusion
Contents
Brownian motion and Diffusion#
import numpy as np
import scipy
from numpy.random import normal, choice, uniform
import ipywidgets as widgets
import matplotlib.pyplot as plt
plt.style.use('fast')
%matplotlib inline
Simulating Brownian motion#
The Brownian motion describes the movement of a particle suspended in a fluid resulting from random collisions with the quick molecules in the fluid (diffusion). More generally, the Brownian motion models a continuous-time random walk, where a particle evolves in space by making independent random steps in all directions.
To simulate brownian motion we take the continuous time limit of random walk and approximate displacements of our particle as normally distributed (binomial->normal, time step->continuous time)
We assume we have started at position \(\mu=0\) and our variance is given by \(\sigma^2=2Dt\) Where D is the diffusion coefficnets which is related to parameters of discree random walk as shown in the lecture.
In the last step we re-wrote brownian motion equation in a convenient way by shifting normally distributed radnom variable by \(\mu\) and scaling it by \(\sigma\)
def brown(T, N, dt=1, D=1):
"""
Creates 3D brownian path given:
time T
N=1 trajecotires
dt=1 timestep
D=1 diffusion coeff
returns np.array with shape (N, T, 3)
"""
nT = int(T/dt) # how many points to sample
dR = np.sqrt(2*D*dt) * np.random.randn(N, nT, 3) # 3D position of brownian particle
R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
return R
Below we proceed to simulate continuous random walk in 1D-3D
We will plot trajecotires and distributions of brownian particle using interactive plotting via ipywidgets and holoviews/plotly interface.
R=brown(T=3000, N=1000)
print(R.shape)
(1000, 3000, 3)
@widgets.interact(t=(10,3000-1))
def brownian_plot(t=10):
fig, ax = plt.subplots(ncols=2)
ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
ax[1].hist(R[:, 10, 0], density=True, color='red')
ax[1].hist(R[:, t, 0], density=True)
ax[1].set_ylim([0,0.1])
ax[0].set_ylim([-200, 200])
ax[0].set_xlim([-200, 200])
fig.tight_layout()
import holoviews as hv
hv.extension('plotly')
plots = []
for i in range(10):
plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
plots.append(plot)
hv.Overlay(plots)
rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]
xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])
hv.Overlay(rw_curve) << ydist << xdist
Diffusion Equation#
The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\)
Diffusion equation:
Formulated empirically as Fick’s laws
This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour
This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.
Can solve numerically by writing derivatives as finite differences!
Can also simulate via random walk!
Diffusion coefficient \(D\), Units \([L^2]/[T]\)
Important special case solution (here written in 1d):
where \(\sigma(t)=\sqrt{2{D}t}\)
density remains Gaussian
Gaussian becomes wider with time
check that this is indeed a solution by plugging into the diffusion equation!
def sigma(t, D = 1):
return np.sqrt(2*D*t)
def gaussian(x, t):
return 1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #
@widgets.interact(t=(1,100,1))
def diffusion(t=0.001):
R = brown(T=101, N=1000)
x = np.linspace(-20, 20, 100)
plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
plt.legend()
plt.ylabel('$p(x)$')
plt.xlabel('$x$')
plt.xlim([-25, 25])
References#
The mighty little books
More in depth
On the applied side